/*
 * case.cpp
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

// uwaga: uzywamy notacji takiej jak w yyt05

#include "case_1D1t_SGA_v2.h"
// MY HEADERS
#include"mros-v42.h"

namespace case_1D1t_SGA_v2
{
global glob;

// *******************************************************************************************  //
// *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY ***
// *******************************************************************************************  //

#define Sqrt sqrt
#define Power pow

/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

struct parameters {
	double n;
	double J;
	double t;
	double U;
	double beta;

	unsigned int L;
};

struct variables {
	double barchiAB;
	double barDelta;
	double mu;
	double mAF;
	double d;
    double Energy;
 };

/*
 *********************************************************
 ***************** Small useful functions ****************
 *********************************************************
*/

inline void clear(parameters &p) {
	p.n = p.J = p.t = p.U = p.beta = 0.; p.L = 0;
}

inline void clear(variables &v) {
	v.barchiAB = v.barDelta  = v.mu = v.mAF = v.d = 0;
	v.Energy = 0;
}

inline void addTo(const variables &q, variables &v, double d=1.) {
	v.barchiAB += q.barchiAB*d;
	v.barDelta += q.barDelta*d;
	v.mu += q.mu*d;
	v.mAF += q.mAF*d;
	v.d += q.d*d;
	v.Energy += q.Energy*d;
}

inline void copyTo(variables &v, const variables &q, double d=1.) {
	v.barchiAB = q.barchiAB*d;
	v.barDelta = q.barDelta*d;
	v.mu = q.mu*d;
	v.mAF = q.mAF*d;
	v.d = q.d*d;
	v.Energy = q.Energy*d;
}

inline void copyTo(parameters &v, const parameters &q) {
	v.n = q.n;
	v.J = q.J;
	v.t = q.t;
	v.U = q.U;
	v.beta = q.beta;
	v.L = q.L;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters  p.n  p.J  p.t  p.tPrim  p.U  p.beta  p.L
	variables v; // variables  v.chiAB  v.chiS  v.chiT  v.barDelta  v.simDelta  v.mu  v.mAF  v.d
	double n, d, m, mu; // contractions

	// Abbreviations - dependend on kx, ky
	double ek;
	double etak;

public:
	double E[4];
	double dEdlbarchiAB[4];
	double dEdlbarDelta[4];
	double dEdln[4];
	double dEdlmAF[4];

	double gt, gs, dgtdm, dgtdd, dgtdn, dgsdm, dgsdd, dgsdn;
	double lbarchiAB, lbarDelta, ln, lmAF;

	Func(const parameters &p, const variables &v) {
		set(p, v);
		refresh();
	}

	void set(const parameters &p, const variables &v) {
		copyTo(this->p, p);
		copyTo(this->v, v);
		refresh();
	}

	void set(const parameters &p) {
		copyTo(this->p, p);
		refresh();
	}

	// Gutzwiller
	double fgt();
	double fgs();
	double fgBDelta();
	double fgm();
	double fdgtdm();
	double fdgtdd();
	double fdgtdn();
	double fdgsdm();
	double fdgsdd();
	double fdgsdn();

    // Energy
    double E1(double ek, double etak);
    double E3(double ek, double etak);

    double dE1dl_barchiAB(double ek, double etak);
    double dE1dl_n(double ek, double etak);
    double dE1dl_mAF(double ek, double etak);
    double dE1dl_barDelta(double ek, double etak);
    double dE1dl_mu(double ek, double etak);

    double dE3dl_barchiAB(double ek, double etak);
    double dE3dl_n(double ek, double etak);
    double dE3dl_mAF(double ek, double etak);
    double dE3dl_barDelta(double ek, double etak);
    double dE3dl_mu(double ek, double etak);

	// DEBUG
	void printLambda() {
		printf("### %25.15e %25.15e %25.15e %25.15e (l)\n", lbarchiAB, lbarDelta, ln, lmAF);
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e (g)\n", gt, gs, dgtdm, dgtdd, dgtdn, dgsdm, dgsdd, dgsdn);
	}

	// REFRESH
	void refresh() {
		n = p.n;
		d = v.d;
		mu = v.mu;
		m = v.mAF;

		// obliczam funkcje Gutzwillara i poszczegolne pochodne
		gt = fgt();
		gs = fgs();
		dgtdn = fdgtdn();
		dgtdd = fdgtdd();
		dgtdm = fdgtdm();
		dgsdn = fdgsdn();
		dgsdm = fdgsdm();
		dgsdd = fdgsdd();

		// obliczanie poszczegolnych lambda:
		lbarchiAB = gt*p.t + (3./8.)*gs*p.J*v.barchiAB;
		lbarDelta = (3./8.)*gs*p.J*v.barDelta;
		ln = 4.*p.t*v.barchiAB*dgtdn \
			+ p.J*(2.*v.mAF*v.mAF + (3./4.)*v.barchiAB*v.barchiAB + (3./4.)*v.barDelta*v.barDelta)*dgsdn;
		lmAF = 4.*p.t*v.barchiAB*dgtdm + 4.*gs*p.J*v.mAF\
			+ p.J*(2.*v.mAF*v.mAF + (3./4.)*v.barchiAB*v.barchiAB + (3./4.)*v.barDelta*v.barDelta)*dgsdm;

		// zeby pozniej mozna bylo wypisac
		glob.lbarchiAB = lbarchiAB;
		glob.lbarDelta = lbarDelta;
		glob.ln = ln;
		glob.lmAF = lmAF;
	}

	//REFRESH XY
	void refresh(double kx, double ky) {
		ek = 2.*(cos(kx) + cos(ky));
		etak = 2.*(cos(kx) - cos(ky));

		refresh();

		E[0]          = E1(ek, etak);
		dEdlbarchiAB[0]  = dE1dl_barchiAB(ek, etak);
		dEdlbarDelta[0] = dE1dl_barDelta(ek, etak);
		dEdln[0]      = dE1dl_n(ek, etak);
		dEdlmAF[0]    = dE1dl_mAF(ek, etak);

		E[1]          = -E[0];
		dEdlbarchiAB[1]  = -dEdlbarchiAB[0];
		dEdlbarDelta[1] = -dEdlbarDelta[0];
		dEdln[1]      = -dEdln[0];
		dEdlmAF[1]    = -dEdlmAF[0];

		E[2]          = E3(ek, etak);
		dEdlbarchiAB[2]  = dE3dl_barchiAB(ek, etak);
		dEdlbarDelta[2] = dE3dl_barDelta(ek, etak);
		dEdln[2]      = dE3dl_n(ek, etak);
		dEdlmAF[2]    = dE3dl_mAF(ek, etak);

		E[3]          = -E[2];
		dEdlbarchiAB[3]  = -dEdlbarchiAB[2];
		dEdlbarDelta[3] = -dEdlbarDelta[2];
		dEdln[3]      = -dEdln[2];
		dEdlmAF[3]    = -dEdlmAF[2];
	}
};

double Func::fgt() {
	return ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
		     (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/
		   (-2*(-m + n/2.)*(m + n/2.) + n);
}

double Func::fgs() {
	return Power(-2*d + n,2)/Power(-2*(-m + n/2.)*(m + n/2.) + n,2);
}

double Func::fgBDelta() {
    //n = 1.-p.sdelta;
    double r = n/2. + v.mAF;
    double w = n/2. - v.mAF;
	double part1 = (n-2.*v.d)/(2.*(n-2.*r*w));
	double part2 = sqrt((1.-w)*(1.-n+v.d)/(1.-r)) + sqrt(w*v.d/r);
	double part3 = sqrt((1.-r)*(1.-n+v.d)/(1.-w)) + sqrt(r*v.d/w);
	return part1*(part2*part2 + part3*part3);
}

double Func::fgm() {
   double r = n/2. + v.mAF;
   double w = n/2. - v.mAF;
   return (n-2.*v.d)/(n-2.*w*r);
}

double Func::fdgsdm() {
	return (-2*(-2*(-m + n/2.) + 2*(m + n/2.))*Power(-2*d + n,2))/Power(-2*(-m + n/2.)*(m + n/2.) + n,3);
}

double Func::fdgsdd() {
	return (-4*(-2*d + n))/Power(-2*(-m + n/2.)*(m + n/2.) + n,2);
}

double Func::fdgsdn() {
	return (-2*(1 - n)*Power(-2*d + n,2))/Power(-2*(-m + n/2.)*(m + n/2.) + n,3) + (2*(-2*d + n))/Power(-2*(-m + n/2.)*(m + n/2.) + n,2);
}

double Func::fdgtdm() {
	return -(((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*(-2*(-m + n/2.) + 2*(m + n/2.))*
	        (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/
	      Power(-2*(-m + n/2.)*(m + n/2.) + n,2)) + ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) +
	        Sqrt((d*(-m + n/2.))/(m + n/2.)))*((-(((1 + d - n)*(1 - m - n/2.))/Power(1 + m - n/2.,2)) - (1 + d - n)/(1 + m - n/2.))/
	         (2.*Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.))) +
	        (d/(-m + n/2.) + (d*(m + n/2.))/Power(-m + n/2.,2))/(2.*Sqrt((d*(m + n/2.))/(-m + n/2.))))*(-2*d + n))/
	    (-2*(-m + n/2.)*(m + n/2.) + n) + ((((1 + d - n)/(1 - m - n/2.) + ((1 + d - n)*(1 + m - n/2.))/Power(1 - m - n/2.,2))/
	         (2.*Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.))) +
	        (-((d*(-m + n/2.))/Power(m + n/2.,2)) - d/(m + n/2.))/(2.*Sqrt((d*(-m + n/2.))/(m + n/2.))))*
	      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/
	    (-2*(-m + n/2.)*(m + n/2.) + n);
}

double Func::fdgtdd() {
	return (-2*(Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
		      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.))))/(-2*(-m + n/2.)*(m + n/2.) + n) +
		   ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
		      ((1 - m - n/2.)/(2.*Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.))*(1 + m - n/2.)) +
		        (m + n/2.)/(2.*(-m + n/2.)*Sqrt((d*(m + n/2.))/(-m + n/2.))))*(-2*d + n))/(-2*(-m + n/2.)*(m + n/2.) + n) +
		   (((1 + m - n/2.)/(2.*(1 - m - n/2.)*Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.))) +
		        (-m + n/2.)/(2.*Sqrt((d*(-m + n/2.))/(m + n/2.))*(m + n/2.)))*
		      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/
		    (-2*(-m + n/2.)*(m + n/2.) + n);
}

double Func::fdgtdn() {
	return -(((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
	        (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(1 - n)*(-2*d + n))/
	      Power(-2*(-m + n/2.)*(m + n/2.) + n,2)) + ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) +
	        Sqrt((d*(-m + n/2.))/(m + n/2.)))*(Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.))))
	     /(-2*(-m + n/2.)*(m + n/2.) + n) + ((Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.)) + Sqrt((d*(-m + n/2.))/(m + n/2.)))*
	      ((((1 + d - n)*(1 - m - n/2.))/(2.*Power(1 + m - n/2.,2)) - (1 + d - n)/(2.*(1 + m - n/2.)) -
	           (1 - m - n/2.)/(1 + m - n/2.))/(2.*Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.))) +
	        (d/(2.*(-m + n/2.)) - (d*(m + n/2.))/(2.*Power(-m + n/2.,2)))/(2.*Sqrt((d*(m + n/2.))/(-m + n/2.))))*(-2*d + n))/
	    (-2*(-m + n/2.)*(m + n/2.) + n) + (((-(1 + d - n)/(2.*(1 - m - n/2.)) +
	           ((1 + d - n)*(1 + m - n/2.))/(2.*Power(1 - m - n/2.,2)) - (1 + m - n/2.)/(1 - m - n/2.))/
	         (2.*Sqrt(((1 + d - n)*(1 + m - n/2.))/(1 - m - n/2.))) +
	        (-(d*(-m + n/2.))/(2.*Power(m + n/2.,2)) + d/(2.*(m + n/2.)))/(2.*Sqrt((d*(-m + n/2.))/(m + n/2.))))*
	      (Sqrt(((1 + d - n)*(1 - m - n/2.))/(1 + m - n/2.)) + Sqrt((d*(m + n/2.))/(-m + n/2.)))*(-2*d + n))/
	    (-2*(-m + n/2.)*(m + n/2.) + n);
}


double Func::E1(double ek, double etak) {
	return -Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lbarDelta,2) + Power(lmAF,2) + 4*Power(ln,2) + 8*ln*mu +
		      4*Power(mu,2) - 4*Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2)))/2.;
}

double Func::E3(double ek, double etak) {
	return -Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lbarDelta,2) + Power(lmAF,2) + 4*Power(ln,2) + 8*ln*mu +
		      4*Power(mu,2) + 4*Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2)))/2.;
}

double Func::dE1dl_barchiAB(double ek, double etak) {
	return -(8*Power(ek,2)*lbarchiAB - (16*Power(ek,2)*lbarchiAB*Power(ln + mu,2))/
		       Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2)))/
		   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lbarDelta,2) + Power(lmAF,2) + 4*Power(ln,2) + 8*ln*mu +
		       4*Power(mu,2) - 4*Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2))));
}

double Func::dE3dl_barchiAB(double ek, double etak) {
	return -(8*Power(ek,2)*lbarchiAB + (16*Power(ek,2)*lbarchiAB*Power(ln + mu,2))/
		       Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2)))/
		   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lbarDelta,2) + Power(lmAF,2) + 4*Power(ln,2) + 8*ln*mu +
		       4*Power(mu,2) + 4*Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2))));
}


double Func::dE1dl_n(double ek, double etak) {
	return -(8*ln + 8*mu - (4*(4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*(ln + mu))/
		       Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2)))/
		   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lbarDelta,2) + Power(lmAF,2) + 4*Power(ln,2) + 8*ln*mu +
		       4*Power(mu,2) - 4*Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2))));
}

double Func::dE3dl_n(double ek, double etak) {
	return -(8*ln + 8*mu + (4*(4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*(ln + mu))/
		       Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2)))/
		   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lbarDelta,2) + Power(lmAF,2) + 4*Power(ln,2) + 8*ln*mu +
		       4*Power(mu,2) + 4*Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2))));
}

double Func::dE1dl_mAF(double ek, double etak) {
	return -(2*lmAF - (4*lmAF*Power(ln + mu,2))/Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2)))/
			   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lbarDelta,2) + Power(lmAF,2) + 4*Power(ln,2) + 8*ln*mu +
			       4*Power(mu,2) - 4*Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2))));
}

double Func::dE3dl_mAF(double ek, double etak) {
	return -(2*lmAF + (4*lmAF*Power(ln + mu,2))/Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2)))/
			   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lbarDelta,2) + Power(lmAF,2) + 4*Power(ln,2) + 8*ln*mu +
			       4*Power(mu,2) + 4*Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2))));
}

double Func::dE1dl_barDelta(double ek, double etak) {
	return (-2*Power(etak,2)*lbarDelta)/
			   Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lbarDelta,2) + Power(lmAF,2) + 4*Power(ln,2) + 8*ln*mu +
			     4*Power(mu,2) - 4*Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2)));
}

double Func::dE3dl_barDelta(double ek, double etak) {
	return (-2*Power(etak,2)*lbarDelta)/
			   Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lbarDelta,2) + Power(lmAF,2) + 4*Power(ln,2) + 8*ln*mu +
			     4*Power(mu,2) + 4*Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2)));
}

double Func::dE1dl_mu(double ek, double etak) {
	return -(8*ln + 8*mu - (4*(4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*(ln + mu))/
		       Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2)))/
		   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lbarDelta,2) + Power(lmAF,2) + 4*Power(ln,2) + 8*ln*mu +
		       4*Power(mu,2) - 4*Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2))));
}

double Func::dE3dl_mu(double ek, double etak) {
	return -(8*ln + 8*mu + (4*(4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*(ln + mu))/
		       Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2)))/
		   (4.*Sqrt(4*Power(ek,2)*Power(lbarchiAB,2) + 4*Power(etak,2)*Power(lbarDelta,2) + Power(lmAF,2) + 4*Power(ln,2) + 8*ln*mu +
		       4*Power(mu,2) + 4*Sqrt((4*Power(ek,2)*Power(lbarchiAB,2) + Power(lmAF,2))*Power(ln + mu,2))));
}

/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	void onePoint(variables &y, int i, int j);

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->set(p);
	}

	void run(variables &w);
};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	f->refresh(x, y);
	if(DEEP_DEBUG) { f->printLambda(); }

	double Etan[4] = { (1-tanh(p.beta*f->E[0]/2.))/2., (1-tanh(p.beta*f->E[1]/2.))/2., (1-tanh(p.beta*f->E[2]/2.))/2., (1-tanh(p.beta*f->E[3]/2.))/2. };

	for(int i=0; i<4; i++) {
		w.barchiAB += Etan[i]*f->dEdlbarchiAB[i];
		w.barDelta += Etan[i]*f->dEdlbarDelta[i];
		w.mu += Etan[i]*f->dEdln[i];
		w.mAF += Etan[i]*f->dEdlmAF[i];
		// w.d w tym miejscu jeszcze nie

        if(f->E[i]>0) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
        else if(f->E[i]<0) { w.Energy += f->E[i]; }
	}
}

void Equa::run(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation
	unsigned int max = p.L/2; // attention: L should be even!

	// SERIAL CASE
	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
	// * edges (2 times)
	for(unsigned int i=1; i<max; i++) {
		onePoint(temp, i, 0);
		onePoint(temp, 0, i);
		onePoint(temp, i, max-i);
	}
   addTo(temp, w, 2);
   clear(temp);
   // * interior area (4 times) - SERIAL CASE
   for(unsigned int i=1; i<max; i++) {
		for(unsigned int j=1; j<i; j++) {
			onePoint(temp, max-i, j);
     	}
	}
	addTo(temp, w, 4);
	clear(temp);
	// * missing corners
	onePoint(w, max, 0);
	//onePoint(w, 0, max);
	onePoint(w, 0, 0);

	// what was out of sumation
	double L2 = (double)(p.L*p.L);
	w.barchiAB /= L2;
	w.barDelta /= L2;
	w.mu /= L2;
	w.mAF /= L2;
    w.Energy /= L2;

	w.barchiAB += 4.*v.barchiAB;
	w.barDelta += 4.*v.barDelta;
	w.mu += -(1-p.n);
	w.mAF += v.mAF;

	// calculate dd
	w.d = p.U - 4.*p.t*v.barchiAB*f->dgtdd \
			+ p.J*(-2.*v.mAF*v.mAF - (3./4.)*v.barchiAB*v.barchiAB - (3./4.)*v.barDelta*v.barDelta)*f->dgsdd;

    v.Energy = w.Energy \
    		- 4.*f->gt*p.t*v.barchiAB + p.U*v.d\
    		+ f->gs*p.J*(-2.*v.mAF*v.mAF - (3./4.)*v.barchiAB*v.barchiAB - (3./4.)*v.barDelta*v.barDelta) \
    		- f->ln*(1-p.n) + f->lmAF*v.mAF - v.mu \
    		+ 4.*(f->lbarchiAB*v.barchiAB + f->lbarDelta*v.barDelta);

    // results
    case_1D1t_SGA_v2::glob.Energy = v.Energy + v.mu*p.n;
	glob.gt = f->fgt();
	glob.gs = f->fgs();
	glob.gBDelta = f->fgBDelta();
	glob.gm = f->fgm();
}

// ***********************************************************************************  //
// *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY ***
// ***********************************************************************************  //


int function(const gsl_vector * x, void *params, gsl_vector * f) {

	struct parameters p;
	p.n = ((Par *) params)->n();
	p.J = ((Par *) params)->J();
	p.t = ((Par *) params)->t();
	p.U = ((Par *) params)->U();
	p.beta = ((Par *) params)->beta();
	p.L = ((Par *) params)->L();

	variables v;
	v.barchiAB = gsl_vector_get (x, 0);
	v.barDelta = gsl_vector_get (x, 1);
	v.mu = gsl_vector_get (x, 2);
	v.mAF = gsl_vector_get (x, 3);
	v.d = gsl_vector_get (x, 4);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4));
	}

    // It's start of the calculations

    Equa E(p, v);
    variables w;
    E.run(w);

    // It's end of the calculations

	gsl_vector_set (f, 0, w.barchiAB);
	gsl_vector_set (f, 1, w.barDelta);
	gsl_vector_set (f, 2, w.mu);
	gsl_vector_set (f, 3, w.mAF);
	gsl_vector_set (f, 4, w.d);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4));
		printf("###  %24.15e# %25.15e# %24.15e# %24.15e# %24.15e | E = %24.15e\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2),
			gsl_vector_get (f, 3), gsl_vector_get (f, 4), case_1D1t_SGA_v2::glob.Energy);
		printf("\n");
	}

    return GSL_SUCCESS;
}

// ************* CLASSES ************* //

// Par

void Par::print(FILE *file) const {
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.6e", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
}

void Par::print_verboten(FILE *file) const {
	fprintf(file, "PARAMETERS: n=%f, J=%f, t=%f, U=%f, beta=%f, L=%d\n", n(), J(), t(), U(), beta(), L());
}

void Par::print_label(FILE *file) const {
	fprintf(file, "n   J   t   U   beta   L\n");
}

// Var

void Var::print(FILE *file) const{
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.12e", tabD[i]); }
}

void Var::print_verboten(FILE *file) const{
	fprintf(file, "VARIABLES:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.12e", tabD[i]); }
	fprintf(file, "\n");
}

void Var::print_label(FILE *file) const {
	fprintf(file, "    barchiAB             barDelta                mu                     mAF                    d");
}

// Res

void Res::set(gsl_multiroot_fsolver *s) {
	for(int i=0; i<n; i++) {
		(*V)[i] = gsl_vector_get (s->x, i);
	}
	this->glob = case_1D1t_SGA_v2::glob;
}

void Res::print(FILE *file) const { V->print(file); fprintf(file, "%21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f\n", glob.Energy, glob.gt, glob.gs, glob.gBDelta, glob.gm, glob.lbarchiAB, glob.lbarDelta, glob.ln, glob.lmAF); }
void Res::print_verboten(FILE *file) const { fprintf(file, "Result: "); V->print_verboten(file); fprintf(file, "%21.15f %21.15f %21.15f %21.15f %21.15f\n", glob.Energy, glob.gt, glob.gs, glob.gBDelta, glob.gm); }
void Res::print_label(FILE *file) const { V->print_label(file); fprintf(file, "                      E                      gt                    gs             gBDelta               gm                 lbarchiAB             lbarDelta             ln                  lmAF\n");}

void example() {

	printf("\nEXAMPLE 1D1t_v2\n");

    Par par(N_PAR);
    par.set_n(0.6);
    par.set_J(1./3.);
    par.set_t(1.);
    par.set_U(12);
    par.set_beta(500.);
    par.set_L(256);
    par.print_verboten();

    Var initV(N_VAR);
    initV.set_barchiAB(3.593609375214e-01);
    initV.set_barDelta(5.102236912890e-04);
    initV.set_mu(1.505490218941e-15);
    initV.set_mAF(-9.029144426464e-01);
    initV.set_d(1.122392620654e-02);

	/*
		GA
		sdelta             BDelta                chi                   m                     d                     mu                      E                      gt                    gs             gBDelta               gm                    BDelta_sc             m_AF                  BDelta_d
		           5.102236912890e-04    3.593609375214e-01    1.505490218941e-15    1.122392620654e-02   -9.029144426464e-01   -1.029980725110211     0.767764116249644     1.890966457949343     0.767764116249644     1.375124160921240     0.000391731441432     0.000000000000002     0.000120601985785

		SGA 1D1t_v1
		sdelta           BDelta                chi                   m                     mu                    d                     lBDelta               lchi                  lm                    ln                      E                      gt                    gs             gBDelta               gm
		          -6.924471251818e-13    3.593602178999e-01   -1.256790616034e-13    4.452148266829e-01    1.122388991608e-02    4.484767526434e-16   -4.795763648411e-15    3.828097687905e-13   -1.348310680017e+00   -1.029982378303645     0.767763802616342     1.890966933224481     0.767763802616342     1.375124333732947


		SGA 1D1t_v2
		sdelta          barchiAB             barDelta                mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm                 lbarchiAB             lbarDelta             ln                  lmAF
		           3.593602179003e-01   -2.818628663568e-13    4.452148266826e-01   -3.646948119761e-13    1.122388991610e-02   -1.029982378303640     0.767763802616518     1.890966933224215     0.767763802616518     1.375124333732850     0.852706088762232    -0.000000000000067    -1.348310680017772     0.000000000001117
	*/
    printf("Starting point: "); initV.print(); printf("\n\n");
    Res res(N_VAR); res.print_label();
    RootGSL root(par, function, res);
    if(root.Solve(initV) == GSL_SUCCESS) { root.print(); }
    printf("Koniec\n");
}

void example_1D_sdelta_run(Par &par, Var &initV) {
	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	initV.shake();
	if(root.Solve(initV) != GSL_SUCCESS) { printf("#D# First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#n           "); res.print_label();
	fflush(stdout);

	LOOPS = 80;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_n(0.996 + (double)i*0.00005);
		printf("%f ", par.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 346;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_n(0.996 - (double)i*0.001);
		printf("%f ", par.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
}

void example_1D_U_run(Par &par, Var &initV) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	int LOOPS;

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      "); res.print_label();
	fflush(stdout);

	initV.shake();
	if(root->Solve(initV) != GSL_SUCCESS) { printf("# First shot is not good.\n"); return; }
	Var old(*res.get_V());

	LOOPS = 100;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_U(10 - (double)i*0.1);
		printf("%f ", par.U());
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); continue; }
		root->print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 100;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_U(10 + (double)i*0.1);
		printf("%f ", par.U());
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); continue; }
		root->print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
	delete root;
}

void example_1D() {
    Par par(N_PAR);
    par.set_n(0.996);
    par.set_J(1./3.);
    par.set_t(1.);
    par.set_U(11.2);
    par.set_beta(500.);
    par.set_L(256);
    par.print_verboten();

	// D m and BDelta
	Var initV_C(N_VAR);
    initV_C.set_barchiAB(1.746895851096951e-01*2);
     initV_C.set_barDelta(1.634351738398281e-01*2);
    initV_C.set_mu(3.787755335751999e+00);
    initV_C.set_mAF(2.309594e-7/2.);
    initV_C.set_d(4.087994655435035e-02*4.087994655435035e-02);

    Var initV_E(N_VAR);
    initV_E.set_barchiAB(1.745446e-01*2);
    initV_E.set_barDelta(1.633045e-01*2);
    initV_E.set_mu(3.777007e+00);
    initV_E.set_mAF(4.793657e-02/2.);
    initV_E.set_d(4.074308e-02*4.074308e-02);

//    example_1D_sdelta_run(par, initV_C);
    example_1D_sdelta_run(par, initV_E);
}

void example_inf() {
    Par par(N_PAR);
    par.set_n(1.-0.4);
    par.set_J(1./3.);
    par.set_t(1.);
    par.set_U(12.);
    par.set_beta(500.);
    par.set_L(256);
    par.print_verboten();

	// I'm choosing start point
	Var initV(N_VAR);

	FILE * pFile;
	pFile = fopen(OUTPUT_FILE, "w");

	// INFINITE
	Res res(N_VAR);
	RootGSL root(par, function, res);
	printf("sdelta      "); res.print_label();
	while(1) {
	    initV.set_barchiAB(frand());
	    initV.set_barDelta(frand());
	    initV.set_mu(16.*frand()-8.);
	    initV.set_mAF(frand());
	    initV.set_d(frand());

		if(root.Solve(initV) == GSL_SUCCESS) {
			printf("       "); root.print(); root.print(pFile); fflush(pFile); fflush(stdout);
		}
	}
	fclose (pFile);
}

void test_energy() {

	struct parameters p;
	p.n = 0.6;
	p.J = 1./3.;
	p.t = 1.;
	p.U = 12.;
	p.beta = 500;
	p.L = 256;

/*
	GA
	sdelta             BDelta                chi                   m                     d                     mu                      E                      gt                    gs             gBDelta               gm                    BDelta_sc             m_AF                  BDelta_d
	           5.102236912890e-04    3.593609375214e-01    1.505490218941e-15    1.122392620654e-02   -9.029144426464e-01   -1.029980725110211     0.767764116249644     1.890966457949343     0.767764116249644     1.375124160921240     0.000391731441432     0.000000000000002     0.000120601985785

	SGA 1D1t_v1
	sdelta           BDelta                chi                   m                     mu                    d                     lBDelta               lchi                  lm                    ln                      E                      gt                    gs             gBDelta               gm
	          -6.924471251818e-13    3.593602178999e-01   -1.256790616034e-13    4.452148266829e-01    1.122388991608e-02    4.484767526434e-16   -4.795763648411e-15    3.828097687905e-13   -1.348310680017e+00   -1.029982378303645     0.767763802616342     1.890966933224481     0.767763802616342     1.375124333732947


	SGA 1D1t_v2
	sdelta          barchiAB             barDelta                mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm                 lbarchiAB             lbarDelta             ln                  lmAF
	           3.593602179003e-01   -2.818628663568e-13    4.452148266826e-01   -3.646948119761e-13    1.122388991610e-02   -1.029982378303640     0.767763802616518     1.890966933224215     0.767763802616518     1.375124333732850     0.852706088762232    -0.000000000000067    -1.348310680017772     0.000000000001117
*/


	struct variables v1;
	v1.barDelta = 5.102236912890e-04;
	v1.barchiAB = 3.593609375214e-01;
	v1.mAF = 1.505490218941e-15;
	v1.mu = 4.452148266829e-01;
	v1.d = 1.122392620654e-02;

//	struct variables v2;
//	v2.barDelta = 2.251021971962e-01;
//	v2.barchiAB = 3.918123843731e-01;
//	v2.mAF = 0;
//	v2.mu = 3.970842610219e+00;
//	v2.d = 2.888476624105e-02;

	struct variables v = v1;

    // It's start of the calculations
    Equa E(p, v);
    variables w;
    E.run(w);
	// It's start of the calculations

    Func *f = new Func(p, v);

	double gt = f->fgt();
	double gs = f->fgs();
	double dgtdn = f->fdgtdn();
	//double dgtdd = f->fdgtdd();
	double dgtdm = f->fdgtdm();
	double dgsdn = f->fdgsdn();
	double dgsdm = f->fdgsdm();
	//double dgsdd = f->fdgsdd();

	double lbarchiAB = gt*p.t + (3./8.)*gs*p.J*v.barchiAB;
	double lbarDelta = (3./8.)*gs*p.J*v.barDelta;
	double ln = 4.*p.t*v.barchiAB*dgtdn \
		+ p.J*(2.*v.mAF*v.mAF + (3./4.)*v.barchiAB*v.barchiAB + (3./4.)*v.barDelta*v.barDelta)*dgsdn;
	double lmAF = 4.*gs*p.J*v.mAF  +4.*p.t*v.barchiAB*dgtdm \
		+ p.J*(2.*v.mAF*v.mAF + (3./4.)*v.barchiAB*v.barchiAB + (3./4.)*v.barDelta*v.barDelta)*dgsdm;

	double W = -4.*f->gt*p.t*v.barchiAB + p.U*v.d\
    		+ f->gs*p.J*(-2.*v.mAF*v.mAF - (3./4.)*v.barchiAB*v.barchiAB - (3./4.)*v.barDelta*v.barDelta);

	double Wtest = -4.*f->gt*p.t*v.barchiAB + p.U*v.d;

    printf("1_v2> %25.15e\n", glob.Energy);
    printf("1_v2> W = %25.15e\n", W);
    printf("1_v2> Wtest = %25.15e\n", Wtest);
    printf("1_v2> gdzie: %25.15e  %25.15e %25.15e %25.15e\n", lbarchiAB- (gt*p.t + (3./8.)*gs*p.J*v.barchiAB), lbarDelta-(3./8.)*gs*p.J*v.barDelta, ln, lmAF-4.*gs*p.J*v.mAF);
    //printf("1_v2> gdzie (test): %25.15e  %25.15e %25.15e  %25.15e\n", dgtdn, dgtdm, dgsdn, dgsdm);

    delete f;
/*
    Par par(N_PAR);
    par.set_n(0.97);
    par.set_J(1./3.);
    par.set_t(1.);
    par.set_U(10.);
    par.set_beta(500.);
    par.set_L(512);
    par.print_verboten();

	//sdelta          barchiAB             barDelta                mu                     mAF                    d                      E
	//                 3.119786e-01         -3.807814e-23          4.999908e+00         -1.624410e-17          4.852439e-01             4.848583	// 0, silne d
	//                -3.918124e-01          2.251022e-01          3.970843e+00         -2.576902e-11          2.888477e-02            -0.368634	// D
	//                -3.859593e-01          1.206862e-01          3.804065e+00         -1.879365e-01          3.084875e-02            -0.363913	// D, m
	//                -4.051486e-01          1.318796e-15          4.143907e+00          7.110104e-13          3.659600e-02            -0.353778	// 0

    double c = 3.859593e-01;
    double Delta = 1.206862e-01;
    double mu = 3.804065e+00;
    double m = 1.879365e-01;
    double d = 3.084875e-02;
   // double zero = 1e-6;
//    double zero = 1e-6;

    Var initV(N_VAR);
    initV.set_barchiAB(c);
    initV.set_barDelta(Delta);
    initV.set_mu(mu);
    initV.set_mAF(m);
    initV.set_d(d);

	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	if(root->Solve(initV) != GSL_SUCCESS) { printf("# First shot is not good.\n"); return; }
	initV.print_label(); printf("\n");
	root->print();
*/
    return;
}

void example_2D_U_run(Par &par, Var &initV, int s1, int s2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	int LOOPS_n, LOOPS_U;
	Var initV_cp(N_VAR); initV_cp = initV;

	double start_U = 11.200000;
	double start_n = 0.996;

	par.set_U(start_U);
	par.set_n(start_n);

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta  U       "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	if(s2>0) { LOOPS_U = 188; } else { LOOPS_U = 112; }
	if(s1>0) { LOOPS_n = 4; } else { LOOPS_n = 546; }

	for(int i=0; i<=LOOPS_n; i++) {
		if(s1>0) { par.set_n(start_n + (double)i*0.001); } else { par.set_n(start_n - (double)i*0.001); }
		par.set_U(start_U);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS_U; j++) {
			if(s2>0) { par.set_U(start_U + (double)j*0.1); } else { par.set_U(start_U - (double)j*0.1); }
			printf("%f %f", par.n(), par.U());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}

		printf("\n");
		fflush(stdout);
	}
	printf("\n#END\n");
	fflush(stdout);
	initV = initV_cp;
}

void example_2D_n_run(Par &par, Var &initV, int s1, int s2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	int LOOPS_n, LOOPS_U;
	Var initV_cp(N_VAR); initV_cp = initV;

	double start_U = 11.200000;
	double start_n = 0.996;

	par.set_U(start_U);
	par.set_n(start_n);

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#U       sdelta        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	if(s2>0) { LOOPS_n = 4; } else { LOOPS_n = 546; }
	if(s1>0) { LOOPS_U = 188; } else { LOOPS_U = 112; }

	for(int i=0; i<=LOOPS_U; i++) {
		if(s1>0) { par.set_U(start_U + (double)i*0.1); } else { par.set_U(start_U - (double)i*0.1); }
		par.set_n(start_n);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS_n; j++) {
			if(s2>0) { par.set_n(start_n + (double)j*0.001); } else { par.set_n(start_n - (double)j*0.001); }
			printf("%f %f", par.U(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}

		printf("\n");
		fflush(stdout);
	}
	printf("\n#END\n");
	fflush(stdout);
	initV = initV_cp;
}

void example_2D() {
	Par par(N_PAR);
	par.set_n(0.996);
	par.set_J(1./3.);
	par.set_t(1.);
	par.set_U(11.2);
	par.set_beta(500);
	par.set_L(256);
	printf("#"); par.print_verboten(); fflush(stdout);

	// I'm choosing start point
	Var initV_A(N_VAR);
	Var initV_C(N_VAR);
	Var initV_D(N_VAR);

	//sdelta          barchiAB             barDelta                mu                     mAF                    d                      E
	//                 3.119786e-01         -3.807814e-23          4.999908e+00         -1.624410e-17          4.852439e-01             4.848583	// 0, silne d
	//                -3.918124e-01          2.251022e-01          3.970843e+00         -2.576902e-11          2.888477e-02            -0.368634	// D
	//                -3.859593e-01          1.206862e-01          3.804065e+00         -1.879365e-01          3.084875e-02            -0.363913	// D, m
	//                -4.051486e-01          1.318796e-15          4.143907e+00          7.110104e-13          3.659600e-02            -0.353778	// 0

	// ************************** //
	// * FOR U=10, sdelta 0.03 * //
	// ************************** //

	// A "zero"
	// 4.051486e-01          4.332663e-13          4.143907e+00          6.121923e-12          3.659600e-02            -0.353775
	initV_A.set_barchiAB(4.051486e-01);
    initV_A.set_barDelta(4.332663e-6);
    initV_A.set_mu(4.143907e+00);
    initV_A.set_mAF(6.121923e-6);
    initV_A.set_d(3.659600e-02);

	// C only BDelta
	// 3.918124e-01         -2.251022e-01          3.970843e+00         -4.071512e-14          2.888477e-02            -0.368634 - tylko D
	// ----> dla n=0.996, U=11.2
	// 11.200000 0.996000          3.493570e-01          3.270549e-01          3.781740e+00         -9.481512e-12          1.647740e-03            -0.244408 - tylko D
	initV_C.set_barchiAB(3.493570e-01);
    initV_C.set_barDelta(3.270549e-01);
    initV_C.set_mu(3.781740e+00);
    initV_C.set_mAF(-9.481512e-7);
    initV_C.set_d(1.647740e-03);
    //11.200000 0.999800    3.393018220983e-01    3.381755872014e-01    3.640576690537e+00   -1.384051603530e-11    7.081316540255e-05   -0.230217228622680     0.001236756226218     3.997267775969175     0.001236756226218     1.999316827311063     0.170771786201380     0.168972297167449    -3.640726341405383    -0.000000000022961
	initV_C.set_barchiAB(3.393018220983e-01);
    initV_C.set_barDelta(3.381755872014e-01);
    initV_C.set_mu(3.640576690537e+00);
    initV_C.set_mAF(1.384051603530e-11);
    initV_C.set_d(7.081316540255e-05);

	// D m and BDelta
    // 3.859593e-01         -1.206862e-01          3.804065e+00          1.879365e-01          3.084875e-02            -0.363913 - koegzystencja! ... tylko, że mniejsza energia!
	initV_D.set_barchiAB(3.859593e-01);
    initV_D.set_barDelta(1.206862e-01);
    initV_D.set_mu(3.804065e+00);
    initV_D.set_mAF(1.879365e-01);
    initV_D.set_d(3.084875e-02);

	// rowniez D m and BDelta
    Var initV_E(N_VAR);
    initV_E.set_barchiAB(1.745446e-01*2);
    initV_E.set_barDelta(1.633045e-01*2);
    initV_E.set_mu(3.777007e+00);
    initV_E.set_mAF(4.793657e-02/2.);
    initV_E.set_d(4.074308e-02*4.074308e-02);
    //10.800000 0.999900 3.359208203853e-01 3.352254918504e-01 3.699063551438e+00 7.432588745386e-02 4.499020940268e-05 -0.230057456173480 0.000701976849694 3.827458095029702 0.000704490338523 1.956389044906381 0.161417334758780 0.160382690305384 -3.699144407719202 0.112144332041226 2 -0.229869859529985 -0.230057456173480
    initV_E.set_barchiAB(3.359208203853e-01);
    initV_E.set_barDelta(3.352254918504e-01);
    initV_E.set_mu(3.699063551438e+00);
    initV_E.set_mAF(7.432588745386e-02);
    initV_E.set_d(4.499020940268e-05);
    //10.800000 0.990000    3.674033265711e-01    2.993470547467e-01    4.175264891393e+00    4.581180681901e-14    7.703597961149e-03   -0.269870075597951     0.095058979843521     3.800084513942092     0.095058979843520     1.949380546210024     0.269579441302706     0.142193013379641    -4.184760466626419     0.000000000000074
//    initV_E.set_barchiAB(3.674033265711e-01);
//    initV_E.set_barDelta(2.993470547467e-01);
//    initV_E.set_mu(4.175264891393e+00);
//    initV_E.set_mAF(4.581180681901e-7);
//    initV_E.set_d(7.703597961149e-03);
    //9.910000 0.999900    3.373820668283e-01    3.343079338179e-01    4.604511128377e+00    7.114447862838e-02    3.458473026590e-04   -0.230154698718505     0.003152635180990     3.836737998991355     0.003153146853718     1.958759300932954     0.164958459678272     0.160331494130408    -4.604592592123128     0.108172811631956
//    initV_E.set_barchiAB(3.373820668283e-01);
//    initV_E.set_barDelta(3.343079338179e-01);
//    initV_E.set_mu(4.604511128377e+00);
//    initV_E.set_mAF(7.114447862838e-02);
//    initV_E.set_d(3.458473026590e-04);
    //11.200000 0.999800    3.362494837648e-01    3.350551731900e-01    3.550356000366e+00    7.339835921625e-02    6.726979905711e-05   -0.230396672929960     0.001202721086514     3.830458471110253     0.001208435395048     1.957155709469804     0.162201431523197     0.160426865804345    -3.550517931441191     0.110931280330479
    initV_E.set_barchiAB(3.362494837648e-01);
    initV_E.set_barDelta(3.350551731900e-01);
    initV_E.set_mu(3.550356000366e+00);
    initV_E.set_mAF(7.339835921625e-02);
    initV_E.set_d(6.726979905711e-05);
    //11.200000 0.999000    3.392625204081e-01    3.333691627054e-01    3.604951676350e+00    6.482744638420e-02    3.526223590791e-04   -0.233273810058959     0.006153339576384     3.855657717913009     0.006174957984428     1.963582877780566     0.169663358977631     0.160669673137426    -3.605771505525063     0.099488867279237
    initV_E.set_barchiAB(3.392625204081e-01);
    initV_E.set_barDelta(3.333691627054e-01);
    initV_E.set_mu(3.604951676350e+00);
    initV_E.set_mAF(6.482744638420e-02);
    initV_E.set_d(3.526223590791e-04);
    //11.200000 0.998000    3.427890550603e-01    3.311595812621e-01    3.667521949055e+00    5.347790335456e-02    7.452306870145e-04   -0.236910581520402     0.012645233549534     3.882813115276369     0.012673588194134     1.970485502427351     0.179018463394719     0.160728845671727    -3.669182896044501     0.083480138157115
    initV_E.set_barchiAB(3.427890550603e-01);
    initV_E.set_barDelta(3.311595812621e-01);
    initV_E.set_mu(3.667521949055e+00);
    initV_E.set_mAF(5.347790335456e-02);
    initV_E.set_d(7.452306870145e-04);
    //11.200000 0.996000    3.490599677709e-01    3.264466118263e-01    3.771796669746e+00    2.669551692120e-02    1.639308549348e-03   -0.244354907311826     0.026518569596700     3.919729730019021     0.026531766116503     1.979830732668584     0.197546161250590     0.159947811204947    -3.775162215550852     0.042729455845163
    initV_E.set_barchiAB(3.490599677709e-01);
    initV_E.set_barDelta(3.264466118263e-01);
    initV_E.set_mu(3.771796669746e+00);
    initV_E.set_mAF(2.669551692120e-02);
    initV_E.set_d(1.639308549348e-03);

//    	printf("# CASE A SGA_v2++\n"); example_2D_U_run(par, initV_A, 1, 1);
//    	printf("# CASE A SGA_v2+-\n"); example_2D_U_run(par, initV_A, 1, -1);
//   	printf("# CASE A SGA_v2-+\n"); example_2D_U_run(par, initV_A, -1, 1);
//    	printf("# CASE A SGA_v2--\n"); example_2D_U_run(par, initV_A, -1, -1);
/*
    	printf("# CASE C SGA_v2++\n"); example_2D_U_run(par, initV_C, 1, 1);
    	printf("# CASE C SGA_v2+-\n"); example_2D_U_run(par, initV_C, 1, -1);
    	printf("# CASE C SGA_v2-+\n"); example_2D_U_run(par, initV_C, -1, 1);
    	printf("# CASE C SGA_v2--\n"); example_2D_U_run(par, initV_C, -1, -1);

    	printf("# CASE D SGA_v2++\n"); example_2D_U_run(par, initV_D, 1, 1);
    	printf("# CASE D SGA_v2+-\n"); example_2D_U_run(par, initV_D, 1, -1);
    	printf("# CASE D SGA_v2-+\n"); example_2D_U_run(par, initV_D, -1, 1);
    	printf("# CASE D SGA_v2--\n"); example_2D_U_run(par, initV_D, -1, -1);
*/

//        	printf("# CASE E SGA_v2++\n"); example_2D_U_run(par, initV_E, 1, 1);
//        	printf("# CASE E SGA_v2+-\n"); example_2D_U_run(par, initV_E, 1, -1);
        	printf("# CASE E SGA_v2-+\n"); example_2D_U_run(par, initV_E, -1, 1);
    		printf("# CASE E SGA_v2--\n"); example_2D_U_run(par, initV_E, -1, -1);

//    	printf("# CASE A SGA_v2++\n"); example_2D_n_run(par, initV_A, 1, 1);
//    	printf("# CASE A SGA_v2+-\n"); example_2D_n_run(par, initV_A, 1, -1);
//    	printf("# CASE A SGA_v2-+\n"); example_2D_n_run(par, initV_A, -1, 1);
//    	printf("# CASE A SGA_v2--\n"); example_2D_n_run(par, initV_A, -1, -1);

    //	printf("# CASE C SGA_v2++\n"); example_2D_n_run(par, initV_C, 1, 1);
    //	printf("# CASE C SGA_v2+-\n"); example_2D_n_run(par, initV_C, 1, -1);
    //	printf("# CASE C SGA_v2-+\n"); example_2D_n_run(par, initV_C, -1, 1);
    //	printf("# CASE C SGA_v2--\n"); example_2D_n_run(par, initV_C, -1, -1);

//    	printf("# CASE D SGA_v2++\n"); example_2D_n_run(par, initV_D, 1, 1);
//    	printf("# CASE D SGA_v2+-\n"); example_2D_n_run(par, initV_D, 1, -1);
//    	printf("# CASE D SGA_v2-+\n"); example_2D_n_run(par, initV_D, -1, 1);
//    	printf("# CASE D SGA_v2--\n"); example_2D_n_run(par, initV_D, -1, -1);


//    	printf("# CASE E SGA_v2++\n"); example_2D_n_run(par, initV_E, 1, 1);
   // 	printf("# CASE E SGA_v2+-\n"); example_2D_n_run(par, initV_E, 1, -1);
//    	printf("# CASE E SGA_v2-+\n"); example_2D_n_run(par, initV_E, -1, 1);
	//	printf("# CASE E SGA_v2--\n"); example_2D_n_run(par, initV_E, -1, -1);

}

void run_1D_U() {
	Par par(N_PAR);
	par.set_n(0.996);
	par.set_J(1./3.);
	par.set_t(1.);
	par.set_U(11.2);
	par.set_beta(500);
	par.set_L(256);
	printf("#"); par.print_label();
	printf("#"); par.print_verboten();

	Var initV_A(N_VAR);
	Var initV_C(N_VAR);
	Var initV_D(N_VAR);
    Var initV_E(N_VAR);

	// (0.03,10) A "zero"
	// 4.051486e-01          4.332663e-13          4.143907e+00          6.121923e-12          3.659600e-02            -0.353775
	initV_A.set_barchiAB(4.051486e-01);
    initV_A.set_barDelta(4.332663e-6);
    initV_A.set_mu(4.143907e+00);
    initV_A.set_mAF(6.121923e-6);
    initV_A.set_d(3.659600e-02);

	// (0.004,11.2) C only BDelta
	// 3.493570e-01          3.270549e-01          3.781740e+00         -9.481512e-12          1.647740e-03            -0.244408 - tylko D
	initV_C.set_barchiAB(3.493570e-01);
    initV_C.set_barDelta(3.270549e-01);
    initV_C.set_mu(3.781740e+00);
    initV_C.set_mAF(-9.481512e-7);
    initV_C.set_d(1.647740e-03);

	// (0.03,10) hihg m and BDelta
    // 3.859593e-01         -1.206862e-01          3.804065e+00          1.879365e-01          3.084875e-02            -0.363913 - koegzystencja! ... tylko, że mniejsza energia!
	initV_D.set_barchiAB(3.859593e-01);
    initV_D.set_barDelta(1.206862e-01);
    initV_D.set_mu(3.804065e+00);
    initV_D.set_mAF(1.879365e-01);
    initV_D.set_d(3.084875e-02);

	// (0.004,11.2) rowniez D m and BDelta
    // 3.490599677709e-01    3.264466118263e-01    3.771796669746e+00    2.669551692120e-02    1.639308549348e-03   -0.244354907311826     0.026518569596700     3.919729730019021     0.026531766116503     1.979830732668584     0.197546161250590     0.159947811204947    -3.775162215550852     0.042729455845163
    initV_E.set_barchiAB(3.490599677709e-01);
    initV_E.set_barDelta(3.264466118263e-01);
    initV_E.set_mu(3.771796669746e+00);
    initV_E.set_mAF(2.669551692120e-02);
    initV_E.set_d(1.639308549348e-03);

	// I'm choosing start point
	Var initV = initV_E;

	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#U           "); res.print_label();
	fflush(stdout);

	LOOPS = 8;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_U(11.2 + ((double)i)*0.1);
		printf("%f ", par.U());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 0;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_U(11.2 - ((double)i)*0.001*5);
		printf("%f ", par.U());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);

}

void run_1D_n() {
	Par par(N_PAR);
	par.set_n(0.996);
	par.set_J(1./3.);
	par.set_t(1.);
	par.set_U(12);
	par.set_beta(500);
	par.set_L(512);
	printf("#"); par.print_label();
	printf("#"); par.print_verboten();

	Var initV_A(N_VAR);
	Var initV_C(N_VAR);
	Var initV_D(N_VAR);
    Var initV_E(N_VAR);
    //Var initV_U30(N_VAR);

	// (0.03,12) A "zero"
    // 12.000000     4.051257658478e-01    6.221108339372e-18    4.132644340635e+00   -1.236737202003e-12    1.355837296602e-02   -0.306856302287931     0.199543723949797     3.562524982530471     0.199543723949797     1.887465226840079     0.379952556687251     0.000000000000000    -4.149069721629376    -0.000000000002724
    initV_A.set_barchiAB(4.051257658478e-01);
    initV_A.set_barDelta(6.221108339372e-12);
    initV_A.set_mu(4.132644340635e0);
    initV_A.set_mAF(1.236737202003e-12);
    initV_A.set_d(1.355837296602e-2);

	// (0.004,12) C "only BDelta"
	// 12.000000     3.472483636367e-01    3.295091443614e-01    3.503684531265e+00   -2.717761697015e-09    9.819780733128e-04   -0.243341123751202     0.020651661296744     3.952557106734233     0.020651661296744     1.988103897369107     0.192216534733522     0.162800463784934    -3.506858740735653    -0.000000004383634
	initV_C.set_barchiAB(3.472483636367e-1);
    initV_C.set_barDelta(3.295091443614e-1);
    initV_C.set_mu(3.503684531265e0);
    initV_C.set_mAF(2.717761697015e-9);
    initV_C.set_d(9.819780733128e-4);

	// (0.03,10) D "hihg m and BDelta"
    // 12.000000     3.596603703871e-01    1.056007132670e-01    3.390488502327e+00    2.730066689828e-01    1.287869746179e-02   -0.324339694341408     0.207928834701003     2.119301912560822     0.215853155267320     1.455782233907538     0.303207448555206     0.027974974199308    -3.574261621214279     0.540472591964946
	initV_D.set_barchiAB(3.596603703871e-1);
    initV_D.set_barDelta(1.056007132670e-1);
    initV_D.set_mu(3.390488502327e0);
    initV_D.set_mAF(2.730066689828e-1);
    initV_D.set_d(1.287869746179e-2);

    // (0.004,12) E "rowniez D m and BDelta"
    //12.000000     3.467588742289e-01    3.287737017049e-01    3.488995221338e+00    3.284140847216e-02    9.772714196983e-04   -0.243352111648176     0.020595039227696     3.918745725537802     0.020620971382518     1.979582209845755     0.190452521249819     0.161047567278146    -3.492325754241237     0.052027806413465
    initV_E.set_barchiAB(3.467588742289e-1);
    initV_E.set_barDelta(3.287737017049e-1);
    initV_E.set_mu(3.488995221338e0);
    initV_E.set_mAF(3.284140847216e-2);
    initV_E.set_d(9.772714196983e-4);

	// I'm choosing start point
	Var initV = initV_E;

	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      "); res.print_label();
	fflush(stdout);

	int MNOZNIK = 1;

	LOOPS = 4*MNOZNIK;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_n(0.996 + ((double)i)*0.001/MNOZNIK);
		printf("%f ", 1.-par.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 446*MNOZNIK;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_n(0.996 - ((double)i)*0.001/MNOZNIK);
		printf("%f ", 1.-par.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
}

void test_L() {
	Par par(N_PAR);
    par.set_n(0.75);
    par.set_J(1./3.);
    par.set_t(1.);
    par.set_U(5);
    par.set_beta(500.);
    par.set_L(128);
	printf("#"); par.print_label();
	printf("#"); par.print_verboten();

    // #U       sdelta                barchiAB             simchiAB                 chiS                  chiT                  barDelta                simDelta                mu                     mAF                    d                      E
    // A: 11.200000 0.970000          2.011346e-01          1.258033e-13         -9.880082e-03         -2.416783e-18         -3.247383e-23          8.555196e-24          4.073413e+00          1.246972e-16          1.453128e-01            -0.320330
    // C: 11.200000 0.970000          1.903212e-01          8.935506e-24         -2.619856e-03         -3.000922e-15          1.320256e-01          9.037352e-24          3.675475e+00         -2.309594e-13          1.131126e-01            -0.344194
    // D: 11.200000 0.970000          1.859003e-01          2.010305e-23          4.103517e-03         -5.301833e-02          5.342205e-02         -4.274270e-23          3.467201e+00          4.775446e-01          1.290175e-01            -0.332657

	// I'm choosing start point
	// A only BDelta
	Var initV_A(N_VAR);
    initV_A.set_barchiAB(2.011346e-01*2);
    initV_A.set_barDelta(-3.247383e-23*2);
    initV_A.set_mu(4.073413e+00);
    initV_A.set_mAF(1.246972e-16/2.);
    initV_A.set_d(1.453128e-01*1.453128e-01);

    //C: 0.004000  1.746895851096951e-01 -7.276282569498597e-17 -1.039785930899697e-03 3.238161681211501e-11 1.634351738398281e-01 -2.033308552312695e-27 3.787755335751999e+00 -2.352695019666226e-09 4.087994655435035e-02   -0.244361775702300     0.026801086418488     3.941603049634213     0.313454890385501     1.828509469934118
	// C only BDelta
	Var initV_C(N_VAR);
    initV_C.set_barchiAB(1.746895851096951e-01*2);
     initV_C.set_barDelta(1.634351738398281e-01*2);
    initV_C.set_mu(3.787755335751999e+00);
    initV_C.set_mAF(2.309594e-7/2.);
    initV_C.set_d(4.087994655435035e-02*4.087994655435035e-02);

	// D m and BDelta
    Var initV_D(N_VAR);
    initV_D.set_barchiAB(1.859003e-01*2);
    initV_D.set_barDelta(5.342205e-02*2);
    initV_D.set_mu(3.467201e+00);
    initV_D.set_mAF(4.775446e-01/2.);
    initV_D.set_d(1.290175e-01*1.290175e-01);

//#U       sdelta     barchiAB     simchiAB        chiS         chiT       barDelta     simDelta     mu            mAF             d              E	        A	     C	       D
//11.200000 0.996000 1.745446e-01 4.111222e-19 -9.708780e-04 -6.585940e-04 1.633045e-01 4.814928e-29 3.777007e+00 4.793657e-02 -4.074308e-02 -0.244423 3 -0.201021 -0.244421 -0.244423
//11.200000 0.998000 1.713966e-01 7.191993e-22 -4.717541e-04 -7.518700e-04 1.656524e-01 1.749248e-31 3.669687e+00 1.058211e-01 -2.737416e-02 -0.236972 3 -0.190431 -0.236916 -0.236972
	// rowniez D m and BDelta
    Var initV_E(N_VAR);
    initV_E.set_barchiAB(1.745446e-01*2);
    initV_E.set_barDelta(1.633045e-01*2);
    initV_E.set_mu(3.777007e+00);
    initV_E.set_mAF(4.793657e-02/2.);
    initV_E.set_d(4.074308e-02*4.074308e-02);

    // #sdelta           BDelta                chi                   m                     mu                    d                     lBDelta               lchi                  lm                    ln                      E                      gt                    gs             gBDelta               gm
    // 0.250000    -6.209804552634e-11    3.895682564385e-01    9.508768873898e-14    8.769895137717e-01    6.187826606100e-02    3.245434167115e-14   -4.172966489299e-14    8.907221238724e-13   -1.422463579468e+00   -1.114840828564147     0.870521193324696     1.784858765356844     0.870521193324696     1.335986064806382
    Var initV_E2(N_VAR);
    initV_E2.set_barchiAB(3.895682564385e-01*2);
    initV_E2.set_barDelta(6.209804552634e-11);
    initV_E2.set_mu(8.769895137717e-01);
    initV_E2.set_mAF(9.508768873898e-14);
    initV_E2.set_d(6.187826606100e-02);

    Var initV = initV_E2;

	Res res(N_VAR);
	RootGSL root(par, function, res);

	initV.shake();
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      "); res.print_label();
	fflush(stdout);


	int LOOPS1 = 10;
	int LOOPS2 = 10;
	int LOOPS3 = 10;
	int LOOPS4 = 10;
	int LOOPS5 = 10;
	int LOOPS6 = 10;

	initV = old;
	for(int i=0; i<LOOPS1; i++) {
		par.set_L(128 + 16*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<LOOPS2; i++) {
		par.set_L(288 + 32*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<LOOPS3; i++) {
		par.set_L(608 + 64*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<LOOPS4; i++) {
		par.set_L(1248 + 128*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<LOOPS5; i++) {
		par.set_L(2528 + 256*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<=LOOPS6; i++) {
		par.set_L(5088 + 512*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}

	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
}

} // end of namespace case2D2t_SGA
